#' @TODO: 利用maf文件，突变分析 v1.0（目前先这么用着，后面可以用table或者其他来进行改进）
#' 
#' 
#' @param maf_file_dir  ：maf文件所在地址，字符串形式 
#'  
#' @param gene_list ：兴趣基因list，字符串向量
#' @param phenotype: 表型文件，maf文件路径，表型文件中，需要有`sample`，`riskgroup`，`Age`，`Stage`，`Gender`等几列，且列名需要一致
#' 在表型文件中，各列表型内容需统一成二分类变量，以便进行注释。二分类变量名，需要和下面的注释一致，或者需要自己改。
#'  
#' @examples 临床信息的排序与重命名
#'  dplyr::rename(Tumor_Sample_Barcode = sample, Group = riskgroup) 
#' 
#'  anno_col1 <- c("#E41A1C", "#377EB8")
#'  names(anno_col1) <- c(sort(unique(maf_data@clinical.data$Group))[1],sort(unique(maf_data@clinical.data$Group))[2])
#'
#'  anno_col2 <- c("#4DAF4A", "#984EA3")
#'  names(anno_col2) <- c(sort(unique(maf_data@clinical.data$Age))[1],sort(unique(maf_data@clinical.data$Age))[2])
#'
#'  anno_col3 <- c("#FF7F00", "#FFFF33")
#'  names(anno_col3) <- c(sort(unique(maf_data@clinical.data$Stage))[1],sort(unique(maf_data@clinical.data$Stage))[2])
#'
#'  anno_col4 <- c("#A65628", "#F781BF")
#'  names(anno_col4) <- c(sort(unique(maf_data@clinical.data$Gender))[1],sort(unique(maf_data@clinical.data$Gender))[2])
#' 
#' 
#' > head(phe_2_1)
#'            sample status      time riskgroup riskscore  Age      Stage
#' 1 TCGA-2J-AAB1-01      1  2.200000      High 0.9894815 >=60 Stage I&II
#' 2 TCGA-2J-AAB4-01      0 24.300000      High 0.8465443  <60 Stage I&II
#' 3 TCGA-2J-AAB6-01      1  9.766667      High 0.7197085 >=60 Stage I&II
#' 4 TCGA-2J-AAB8-01      0  2.666667      High 0.3918442 >=60 Stage I&II
#' 5 TCGA-2J-AAB9-01      1 20.900000       Low 0.2620438 >=60 Stage I&II
#' 6 TCGA-2J-AABA-01      1 20.233333       Low 0.2359602  <60 Stage I&II
#'   Gender
#' 1   male
#' 2   male
#' 3   male
#' 4   male
#' 5 female
#' 6   male
#' 
#'@param return_file_style: maf文件数据，生成本地图片
#'@param return_file: 
#'
#'@Author *WYK*
######
maf_analize <- function(phenotype = NULL,Feature = NULL, maf_file_dir = NULL, gene_list = NULL, topN = 20, output_dir = './') {
  library(maftools)

  cli <- phenotype %>% #, age, stage, gender
    dplyr::rename(Tumor_Sample_Barcode = sample, Group = riskgroup) %>%
    mutate(Tumor_Sample_Barcode = substr(Tumor_Sample_Barcode, 1, 12))

  maf_data <- read.maf(
    maf = maf_file_dir,
    clinicalData = cli,
    isTCGA = T,
    verbose = T
  )

  vc_cols <- RColorBrewer::brewer.pal(n = 8, name = "Paired")
  names(vc_cols) <- c(
    "Frame_Shift_Del",
    "Missense_Mutation",
    "Nonsense_Mutation",
    "Multi_Hit",
    "Frame_Shift_Ins",
    "In_Frame_Ins",
    "Splice_Site",
    "In_Frame_Del"
  )

  Gene_mut_res <- maftools::getGeneSummary(x = maf_data) %>%
    as.data.frame() %>%
    arrange(desc(total))

  top_mut_gene <- Gene_mut_res %>%
    filter(Hugo_Symbol %in% gene_list) %>%
    arrange(desc(total)) %>%
    dplyr::select(Hugo_Symbol) %>%
    pull() %>%
    .[1:topN]

  if (!dir.exists(sprintf("%s/snv/", output_dir))) {
    dir.create(sprintf("%s/snv/", output_dir), recursive = T)
  }

  pdf(file = sprintf("%s/snv/SNV_top_genelist.pdf",output_dir), width = 7, height = 7.5)
  oncoplot(
    maf = maf_data,
    top = topN,
    fontSize = .8,
    draw_titv = T,
    legendFontSize = 1.5,
    # clinicalFeatures = c("Age", "Stage", "Gender"),
    colors = vc_cols,
    annotationFontSize = 1.2,
    # annotationColor = anno_col,
    anno_height = 2.6,
    genes = top_mut_gene,
    sortByAnnotation = F
  )
  dev.off()

 
  print(head(maf_data@clinical.data))


  anno_col1 <- c("#E41A1C", "#377EB8")
  names(anno_col1) <- c(sort(unique(maf_data@clinical.data$Group))[1],sort(unique(maf_data@clinical.data$Group))[2])

  anno_col2 <- c("#4DAF4A", "#984EA3")
  names(anno_col2) <- c(sort(unique(maf_data@clinical.data$Age))[1],sort(unique(maf_data@clinical.data$Age))[2])

  anno_col3 <- c("#FF7F00", "#FFFF33")
  names(anno_col3) <- c(sort(unique(maf_data@clinical.data$Stage))[1],sort(unique(maf_data@clinical.data$Stage))[2])

  anno_col4 <- c("#A65628", "#F781BF")
  names(anno_col4) <- c(sort(unique(maf_data@clinical.data$Gender))[1],sort(unique(maf_data@clinical.data$Gender))[2])

  anno_col <- list(Group = anno_col1, Age = anno_col2, Stage = anno_col3) # , Gender = anno_col4


  pdf(file = sprintf("%s/snv/SNV_HighLow.pdf",output_dir), width = 7, height = 7.5)
  oncoplot(
    maf = maf_data,
    top = topN,
    fontSize = .8,
    draw_titv = F,
    legendFontSize = 1.5,
    clinicalFeatures = c("Group", "Age", "Stage", "Gender"),
    colors = vc_cols,
    annotationFontSize = 1.2,
    annotationColor = anno_col,
    anno_height = 2.6,
    # genes = inf_gene$Symbol,
    sortByAnnotation = TRUE
  )
  dev.off()

  maf_data_high <- maf_data
  maf_data_low <- maf_data

  maf_data_high <- subsetMaf(maf_data, tsb = {
    maf_data@clinical.data %>%
      filter(Group == "High") %>%
      dplyr::select(Tumor_Sample_Barcode) %>%
      pull()
  })

  maf_data_low <- subsetMaf(maf_data, tsb = {
    maf_data@clinical.data %>%
      filter(Group == "Low") %>%
      dplyr::select(Tumor_Sample_Barcode) %>%
      pull()
  })

# anno_col_1 <- anno_col
# anno_col_1[[1]] <- anno_col[[1]][1]
# anno_col_2 <- anno_col
# anno_col_2[[1]] <- anno_col[[1]][2]

  
# top20_common_Gene <- intersect(
#   maf_data_high@gene.summary %>% .[1:20, ] %>% pull(1),
#   maf_data_low@gene.summary %>% .[1:20, ] %>% pull(1)
# )

# maftools::coOncoplot(
#   m1 = maf_data_high, 
#   m2 = maf_data_low,
#   genes = top20_common_Gene,
#   clinicalFeatures1 = c("Group", "Age", "Stage", "Gender"),
#   clinicalFeatures2 = c("Group", "Age", "Stage", "Gender"),
#   anno_height = 2.6,
#   annotationColor1 = anno_col_1,
#   annotationColor2 = anno_col_2,
#   sortByAnnotation1 = T,
#   sortByAnnotation2 = T,
#   m1Name = 'High Score',
#   m2Name = 'Low Score'
# )

  pdf(file = sprintf("%s/snv/SNV_High_group.pdf",output_dir), width = 5, height = 7.5)
  oncoplot(
    maf = maf_data_high,
    top = topN,
    fontSize = .8,
    draw_titv = F,
    legendFontSize = 1.5,
    clinicalFeatures = c("Group", "Age", "Stage", "Gender"),
    colors = vc_cols,
    annotationFontSize = 1.2,
    annotationColor = anno_col,
    anno_height = 2.6,
    # genes = c('WTAP','RBM15'),
    sortByAnnotation = TRUE
  )
  dev.off()


  pdf(file = sprintf("%s/snv/SNV_Low_group.pdf",output_dir), width = 5, height = 7.5)
  oncoplot(
    maf = maf_data_low,
    top = topN,
    fontSize = .8,
    draw_titv = F,
    legendFontSize = 1.5,
    clinicalFeatures = c("Group", "Age", "Stage", "Gender"),
    colors = vc_cols,
    annotationFontSize = 1.2,
    annotationColor = anno_col,
    anno_height = 2.6,
    # genes = c('WTAP','RBM15'),
    sortByAnnotation = TRUE
  )
  dev.off()

  mut_matrix <- maftools::mutCountMatrix(maf_data) %>% as.data.frame()

  mut_matrix_long <- mut_matrix %>%
    tibble::rownames_to_column("gene") %>%
    tidyr::pivot_longer(cols = -gene, names_to = "sample", values_to = "mut") %>%
    full_join(maf_data@clinical.data %>% select(Tumor_Sample_Barcode, Group) %>% rename(sample = Tumor_Sample_Barcode)) %>%
    mutate(mut = ifelse(mut > 0, "mut", "wt"))

  genelist <- unique(mut_matrix_long %>% pull(gene)) %>% na.omit() %>% as.character()

  chisq.test.p <- map_dbl(genelist, function(x) {
    # x <- "ZNF789"
    mut_matrix_long %>%
      filter(gene == x) %>%
      select(mut, Group) %>%
      table() %>%
      chisq.test(., simulate.p.value = T) %>%
      .$p.val
  })

  chisq.test_res <- tibble(gene = genelist, chisq.test.p = chisq.test.p) %>% arrange(chisq.test.p)

  output_dir <- "/Pub/Users/wangyk/project/Poroject/P220329001_GDT220112001_CHOL_hypoxia/out/2."

  write_tsv(x = chisq.test_res, file = sprintf("%s/snv/chisq.test_res.txt",output_dir))

  tmb_1 <- tmb(maf = maf_data) %>%
    left_join(., maf_data@clinical.data %>% dplyr::select(Tumor_Sample_Barcode, Group))

  p <- ggplot(tmb_1 %>% filter(total_perMB_log != "-Inf" & Group != "NA"), aes(Group, total_perMB_log)) +
    geom_boxplot(aes(group = Group, fill = Group), color = "black", alpha = 0.1, width = 0.4, notch = T) +
    geom_jitter(aes(colour = Group), position = position_jitter(0.12), alpha = 0.8) +
    geom_violin(aes(group = Group, fill = Group), color = "black", alpha = 0.1, width = 0.55) +
    ggpubr::stat_compare_means() +
    egg::theme_article(base_size = 14) +
    theme(legend.position = "none") +
    scale_fill_manual(values = c("High" = "#E41A1C", "Low" = "#377EB8")) +
    scale_color_manual(values = c("High" = "#E41A1C", "Low" = "#377EB8")) +
    labs(x = "Risk Group", y = "total perMB log")

  ggsave(sprintf("%s/snv/High_low_box_p.pdf",output_dir), width = 4, height = 4)

  return(maf_data)
}


